Organising principles among vital rates: Reduced models (without time interval)
Data
Organising all tables of models results for all 21 forest plo in one object per vital rate, quadrat size scale.
See table S1 for forest plot names.
Growth
Import results data to list
path = here("models_outputs/models_notime/all", "grow", "table")
files = list.files(path)
all = lapply(files, function(x) get(assign(x, load(paste0(path, "/", x)))))
tudo <- bind_rows(all)Remove intercept to include again as a separate column
intercepts <- tudo %>% filter(term == "intercept")
names(intercepts)[8] = "intercept"
grow <- tudo %>% filter(term != "intercept") %>%
left_join(intercepts[,c(1:4,8)], by= c("data", "fplot", "time", "q_size")) %>%
mutate(q_size = fct_relevel(q_size,"quad_5", "quad_10","quad_20",
"quad_50", "quad_100")) %>%
group_by(data, fplot, time, q_size) %>%
mutate(VPC = variance/sum(variance)) %>%
as.data.frame()Divide the sd of the terms by the mean growth (intercept) for each dataset.
grow$stdev = grow$Estimate/grow$interceptNumber of growth models per florest plot and quadrat scale.Some forest plots have more than 1 model because we ran one model per census interval.
table(grow$fplot,grow$q_size)/4##
## quad_5 quad_10 quad_20 quad_50 quad_100
## ama 1 1 1 1 1
## bci 6 6 6 6 6
## edo 2 2 2 2 2
## fus 3 3 3 3 3
## idc 1 1 1 1 1
## kor 1 1 1 1 1
## lam 3 3 3 3 3
## ldw 1 1 1 1 1
## len 2 2 2 2 2
## lpl 1 1 1 1 1
## luq 3 3 3 3 3
## mos 2 2 2 2 2
## pas 4 4 4 4 4
## scbi 1 1 1 1 1
## serc 1 1 1 1 1
## sin 2 2 2 2 2
## ucsc 2 2 2 2 2
## wab 2 2 2 2 2
## wfdp 1 1 1 1 1
## wyw 3 3 3 3 3
## zof 1 1 1 1 1
Mortality
Import results data to list
path = here("models_outputs/models_notime/all", "mort", "table")
files = list.files(path)
all = lapply(files, function(x) get(assign(x, load(paste0(path, "/", x)))))
tudo <- bind_rows(all)Remove intercept to include again as a separate column
intercepts <- tudo %>% filter(term == "intercept")
names(intercepts)[8] = "intercept"
mort <- tudo %>% filter(term != "intercept") %>%
left_join(intercepts[,c(1:4,8)], by= c("data", "fplot", "time", "q_size")) %>%
mutate(q_size = fct_relevel(q_size,"quad_5", "quad_10","quad_20",
"quad_50", "quad_100")) %>%
group_by(data, fplot, time, q_size) %>%
mutate(VPC = variance/sum(variance)) %>%
as.data.frame()
mort$stdev <- mort$EstimateNumber of mortality models per florest plot and quadrat scale.Some forest plots have more than 1 modelo because we ran one model per census interval.
table(mort$fplot,mort$q_size)/4##
## quad_5 quad_10 quad_20 quad_50 quad_100
## ama 1 1 1 1 1
## bci 7 7 7 7 7
## edo 2 2 2 2 2
## fus 3 3 3 3 3
## idc 1 1 1 1 1
## kor 1 1 1 1 1
## lam 3 3 3 3 3
## ldw 1 1 1 1 1
## len 2 2 2 2 2
## lpl 1 1 1 1 1
## luq 3 3 3 3 3
## mos 2 2 2 2 2
## pas 5 5 5 5 5
## scbi 1 1 1 1 1
## serc 1 1 1 1 1
## sin 2 2 2 2 2
## ucsc 2 2 2 2 2
## wab 2 2 2 2 2
## wfdp 1 1 1 1 1
## wyw 3 3 3 3 3
## zof 1 1 1 1 1
Recruitment
Import results data to list
path = here("models_outputs/models_notime/all", "rec", "table")
files = list.files(path)
all = lapply(files, function(x) get(assign(x, load(paste0(path, "/", x)))))
tudo <- bind_rows(all)Remove intercept to include again as a separate column
intercepts <- tudo %>% filter(term == "intercept")
names(intercepts)[8] = "intercept"
rec <- tudo %>% filter(term != "intercept") %>%
left_join(intercepts[,c(1:4,8)], by= c("data", "fplot", "time", "q_size")) %>%
mutate(q_size = fct_relevel(q_size,"quad_5", "quad_10","quad_20",
"quad_50", "quad_100")) %>%
group_by(data, fplot, time, q_size) %>%
mutate(VPC = variance/sum(variance)) %>%
as.data.frame()
rec$stdev <- rec$EstimateNumber of recruitment models per florest plot and quadrat scale. Some forest plots have more than 1 modelo because we ran one model per census interval.
table(rec$fplot,rec$q_size)/4##
## quad_5 quad_10 quad_20 quad_50 quad_100
## bci 7 7 7 7 7
## edo 2 2 2 2 2
## fus 3 3 3 3 3
## idc 1 1 1 1 1
## kor 1 1 1 1 1
## lam 3 3 3 3 3
## ldw 1 1 1 1 1
## len 2 2 2 2 2
## lpl 1 1 1 1 1
## luq 3 3 3 3 3
## mos 2 2 2 2 2
## pas 5 5 5 5 5
## scbi 1 1 1 1 1
## serc 1 1 1 1 1
## sin 2 2 2 2 2
## ucsc 2 2 2 2 2
## wab 2 2 2 2 2
## wfdp 1 1 1 1 1
## wyw 3 3 3 3 3
## zof 1 1 1 1 1
Excluding WYW recruitment data because it has very low recruitment rates in all 3 intervals, 0.0001, 0.0004, 0.0007, respectivelly.
rec <- rec %>% filter(fplot!="wyw")Combining models’ results
All vita rates
vital <- bind_rows(list(grow = grow, mort = mort, rec = rec),
.id="vital") %>%
mutate(q_size = fct_recode(q_size,`5`="quad_5",`10`="quad_10",`20`="quad_20",
`50`="quad_50", `100`="quad_100")) %>%
mutate(q_size = fct_relevel(q_size,"5", "10", "20", "50", "100"),
term = fct_relevel(term, "quadrat", "quadrat:sp", "sp", "Residual")) Plot information
Richness by rarefaction at the 6ha plot size (sampling unit 20x20m).
load(here("data", "rarefaction.curves.Rdata"))
rich.rare <- rare[rare$sites ==150,2:3]
colnames(rich.rare)[1] <- "richness.rarefaction"
vital <- vital %>% left_join(rich.rare, by="fplot")Averages
To get ONE value per forest-vital, first I calculate the mean of the posterior variances for the plots >1 interval. Then I use this value as the Variance and calculate the Variance Partition Coefficient (VPC).
mvital <- vital %>% group_by(vital,fplot, term, q_size) %>%
summarise(variance = mean(variance),
stdev = mean(stdev),
richness.rarefaction = mean(richness.rarefaction),
ntrees = mean(ntrees)) %>%
group_by(vital, fplot, q_size) %>%
mutate(VPC = variance/sum(variance))Mean VPC across all plots
avital <- mvital %>% group_by(vital, term, q_size) %>%
summarise(mean.stdev = mean(stdev),
mean.var = mean(variance),
mean.VPC = mean(VPC),
sd.VPC = sd(VPC)) TOTAL VPC across plots
tvital <- mvital %>% group_by(vital, term, q_size) %>%
summarise(variance = sum(variance),
stdev = sqrt(variance)) %>%
group_by(vital, q_size) %>%
mutate(VPC = variance/sum(variance))Saving objects for models results
save(mvital, vital, file= here("models_outputs", "all_results_notime.Rdata"))Standard deviations
Stacked SD for all forest plots
mvital %>%
ggplot(aes(x=fct_reorder(fplot,richness.rarefaction,max), y=stdev, fill=term,col=term)) +
geom_col(alpha=0.8) +
scale_fill_manual(values= cores)+ scale_color_manual(values= cores)+
facet_grid(q_size~vital, scales="free") +
xlab("Number of species -->") +
theme(legend.position="top",
axis.text.x = element_text(angle=45, hjust=1))Standard deviations x species richness
5x5m quadrat scale
overallSD <- mvital %>% group_by(vital, fplot, q_size, term, richness.rarefaction) %>%
summarise(stdev = mean(stdev)) %>%
group_by(vital, fplot, q_size,richness.rarefaction) %>%
summarise(ALLstdev = sum(stdev)) %>%
filter(q_size=="5")Models SD with log(richness) per vital rate
mg <- lm(ALLstdev ~ log(richness.rarefaction),
data=overallSD[overallSD$vital == "grow", ])
mm <- lm(ALLstdev ~ log(richness.rarefaction),
data=overallSD[overallSD$vital == "mort", ])
mr <- lm(ALLstdev ~ log(richness.rarefaction),
data=overallSD[overallSD$vital == "rec", ])
coefs <- bind_rows(list(grow=broom::tidy(mg),
mort=broom::tidy(mm),
rec =broom::tidy(mr)), .id="vital")
vals = c(paste0("beta = ", round(coefs[2,3],2), "; p = ", round(coefs[2,6],3)),
paste0("beta = ", round(coefs[4,3],2), "; p = ", round(coefs[4,6],3)),
paste0("beta = ", round(coefs[6,3],2), "; p = ", round(coefs[6,6],3)))
anota <- data.frame(vital=c("grow", "mort", "rec"),
x=10, y=8,
lab = vals)Overall standard deviations decreased with richness for recruitment.
over <- overallSD %>%
ggplot(aes(x=richness.rarefaction, y=ALLstdev))+
geom_point() +
facet_grid(~vital, labeller = labvit) +
scale_x_log10() +
geom_smooth(method="lm", aes(linetype=vital)) +
scale_linetype_manual(values=c("dashed", "dashed", "solid"))+
theme(legend.position = "none",
panel.background = element_rect(colour="lightgray"))+
geom_text(data=anota, aes(x=x,y=y, label=lab), hjust=0) +
xlab("Rarefied species richness (log)") +
ylab("Standard deviation")
overpng(here("figures", "FIG_S5.3_SDxRicnhess.png"), height = 350, width=900, res=100)
over
dev.off()## quartz_off_screen
## 2
Relationship between standard deviation for each term and species richness.
mgrow <- lm(stdev ~ -1 + log(richness.rarefaction)*term,
data=mvital[mvital$vital == "grow",])
mmort <- lm(stdev ~ -1 + log(richness.rarefaction)*term,
data=mvital[mvital$vital == "mort" ,])
mrec <- lm(stdev ~ -1+log(richness.rarefaction)*term,
data=mvital[mvital$vital == "rec" ,])
broom::glance(mgrow)## # A tibble: 1 × 12
## r.squared adj.r.squ…¹ sigma stati…² p.value df logLik AIC BIC devia…³
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.933 0.932 0.204 722. 5.55e-237 8 75.6 -133. -96.9 17.2
## # … with 2 more variables: df.residual <int>, nobs <int>, and abbreviated
## # variable names ¹​adj.r.squared, ²​statistic, ³​deviance
broom::glance(mmort)## # A tibble: 1 × 12
## r.squared adj.r.squ…¹ sigma stati…² p.value df logLik AIC BIC devia…³
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.946 0.945 0.207 909. 2.12e-256 8 69.0 -120. -83.6 17.7
## # … with 2 more variables: df.residual <int>, nobs <int>, and abbreviated
## # variable names ¹​adj.r.squared, ²​statistic, ³​deviance
broom::glance(mrec)## # A tibble: 1 × 12
## r.squared adj.r.squ…¹ sigma stati…² p.value df logLik AIC BIC devia…³
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.942 0.941 0.269 756. 7.34e-225 8 -36.0 90.0 125. 26.9
## # … with 2 more variables: df.residual <int>, nobs <int>, and abbreviated
## # variable names ¹​adj.r.squared, ²​statistic, ³​deviance
bonferroni correction -> 0.05/3 = 0.0166667
Taking coefficients for figure
coefs <- bind_rows(list(grow=broom::tidy(mgrow),
mort=broom::tidy(mmort),
rec =broom::tidy(mrec)), .id="vital") %>%
filter(! term %in% c("termquadrat", "termquadrat:sp", "termsp",
"termResidual")) %>%
mutate(sig = ifelse(p.value<0.01666667, "yes", "no"),
beta = paste0("beta = ", round(estimate,2), "; p = ", round(p.value,4)),
xis = 150,
y= 3)
coefs$term[coefs$term == "log(richness.rarefaction)"] <- "quadrat"
coefs$term[coefs$term == "log(richness.rarefaction):termsp"] <- "sp"
coefs$term[coefs$term == "log(richness.rarefaction):termquadrat:sp"] <- "quadrat:sp"
coefs$term[coefs$term == "log(richness.rarefaction):termResidual"] <- "Residual"
coefs <- coefs %>% mutate(term = fct_relevel(term, "quadrat","quadrat:sp", "sp", "Residual"))
coefs$beta[coefs$term=="Residual" & coefs$vital !="grow"] <- ""Figure
mvital$line <- "dashed"
mvital$line[mvital$vital %in% c("grow", "rec") & mvital$term == "sp" ] <- "solid"
mvital$line[mvital$vital == "mort" & mvital$term == "quadrat:sp" ] <- "solid"
fig <- mvital %>% filter(q_size=="5") %>%
ggplot(aes(x=richness.rarefaction, y=stdev))+
facet_grid(vital~term, labeller = grpvit) +
geom_point()+
scale_x_log10() +
geom_smooth(method="lm", aes(linetype=line)) +
scale_linetype_manual(values=c("dashed", "solid"))+
geom_text(data=coefs, aes(x=xis,y=y, label=beta))+
theme(legend.position = "none",
panel.background = element_rect(colour="lightgray")) +
xlab("Rarefied species richness (log)") + ylab("Standard deviation")
figpng(here("figures", "FIG_S5.4_SD_regressions.png"), height = 700, width=900, res=100)
fig
dev.off()## quartz_off_screen
## 2
Variance Partition Coefficients VPC
Stacked VPC for all forest plots at the 5x5 m quadrat scale.
vital %>% group_by(vital, fplot, q_size, term) %>%
summarise(VPC = mean(VPC),
richness.rarefaction = mean(richness.rarefaction),
ntrees = mean(ntrees)) %>%
filter(q_size == "5") %>%
ggplot(aes(x=fct_reorder(fplot, richness.rarefaction,max ),
y=VPC, fill=term, col=term)) +
geom_col(alpha=0.8) +
scale_fill_manual(values= cores)+
scale_color_manual(values= cores)+
facet_wrap(~vital,ncol=1, labeller = labvit) +
xlab("- Number of species +") +
theme(axis.text.x = element_text(angle=45, hjust=1),
legend.position="none")ah <- mvital %>% filter(q_size == "5") %>%
group_by(vital,term) %>%
summarise(VPC_mean = as.character(round(mean(VPC), 2))) %>%
mutate(x_val = as.numeric(VPC_mean),
y_val = 0.7:3.7)
ah[5,3] <- c("0.10")Figure 1b: VPCs for all forest plots and vital rates
At the 5x5m quadrat size scale.
vpc.fig1 <- mvital %>% filter(q_size == "5") %>%
ggplot(aes(x=VPC, y=rev(term), col=term)) +
facet_grid(~vital, labeller = labvit) +
geom_jitter(height = 0.1, size=3, alpha=0.6) +
scale_fill_manual(values=cores, labels=grplab) +
scale_color_manual(values=cores, labels=grplab) +
scale_y_discrete(labels = c( "Residual","Species","Species x Space", "Space")) +
ylab("") +
stat_summary(fun=mean, fun.min=mean, fun.max=mean, fatten=1, col="black",
geom="crossbar", width=0.4)+
geom_vline(xintercept=c(0,0.2,0.4,0.6), col='gray', linetype="dashed", alpha=0.6)+
geom_text(data=ah, aes(x=x_val, y=rev(y_val), label=VPC_mean),
hjust=0.5,size=5,
col="black")+
theme(legend.position = "none",
panel.background = element_rect(colour="lightgray"),
axis.text.x = element_text(angle=45, hjust=1, size=16))
vpc.fig1avitsum <- avital %>%
filter(q_size == "5") %>%
select(vital,term,mean.VPC, sd.VPC) %>%
pivot_wider(names_from = vital, values_from = c(mean.VPC, sd.VPC)) %>%
select(term, mean.VPC_grow, sd.VPC_grow, mean.VPC_mort, sd.VPC_mort,
mean.VPC_rec, sd.VPC_rec)
avitsum %>% mutate_if(is.numeric, round, digits=2) %>%
htmlTable(header=c("" ,rep(c("Mean", "SD"), 3)),
cgroup = c("", "Growth", "Mortality", "Recruitment"),
n.cgroup =c(1,2,2,2),
caption = "Average VPC values for models at the 5x5m quadrat grain size.")| Average VPC values for models at the 5x5m quadrat grain size. | ||||||||||
| Â | Growth | Â | Mortality | Â | Recruitment | |||||
|---|---|---|---|---|---|---|---|---|---|---|
| Â | Mean | SD | Â | Mean | SD | Â | Mean | SD | ||
| 1 | quadrat | Â | 0.04 | 0.02 | Â | 0.1 | 0.06 | Â | 0.18 | 0.09 |
| 2 | quadrat:sp | Â | 0.13 | 0.06 | Â | 0.14 | 0.06 | Â | 0.15 | 0.06 |
| 3 | sp | Â | 0.29 | 0.12 | Â | 0.29 | 0.12 | Â | 0.36 | 0.15 |
| 4 | Residual | Â | 0.54 | 0.12 | Â | 0.48 | 0.1 | Â | 0.31 | 0.11 |
Figure 1c: Average VPCs
5x5 m quadrat size scale.
tot.vpc.fig <- tvital %>% filter(q_size == "5") %>%
ggplot(aes(x=vital, fill=term, col=term, y=VPC)) +
geom_col(alpha=0.7,width = 0.6) +
facet_grid(~vital, scales="free_x", labeller = labvit) +
xlab("") + ylab("Mean VPC") +
scale_fill_manual(values=cores, labels=grplab) +
scale_color_manual(values=cores, labels=grplab) +
theme(axis.text.x = element_blank(),
axis.ticks = element_blank())
tot.vpc.figFigure 3: Total VPC per vital rates across spatial grain sizes
barvpc <- tvital %>%
ggplot(aes(x=q_size, fill=term, col=term, y=VPC)) +
geom_col(alpha=0.7) +
facet_grid(~vital, scales="free_x", labeller = labvit) +
xlab("Spatial grain (quadrat size)") + ylab("VPC") +
scale_fill_manual(values=cores, labels=grplab) +
scale_color_manual(values=cores, labels=grplab) +
theme(legend.position = "none") +
labs(tag="(a)")
barvpcpng(here("figures", "FIG_3_VPC_notime_scales.png"), height = 550, width=1080, res=130)
barvpc
dev.off()## quartz_off_screen
## 2
tvital %>% select(vital, term, q_size, VPC) %>%
pivot_wider(names_from = q_size, values_from = VPC) %>%
mutate_if(is.numeric, round, digits=2) %>%
htmlTable()| vital | term | 5 | 10 | 20 | 50 | 100 | |
|---|---|---|---|---|---|---|---|
| 1 | grow | quadrat | 0.02 | 0.02 | 0.03 | 0.02 | 0.02 |
| 2 | grow | quadrat:sp | 0.14 | 0.11 | 0.08 | 0.08 | 0.07 |
| 3 | grow | sp | 0.33 | 0.32 | 0.32 | 0.28 | 0.26 |
| 4 | grow | Residual | 0.51 | 0.54 | 0.57 | 0.62 | 0.65 |
| 5 | mort | quadrat | 0.1 | 0.09 | 0.06 | 0.04 | 0.03 |
| 6 | mort | quadrat:sp | 0.15 | 0.12 | 0.11 | 0.08 | 0.05 |
| 7 | mort | sp | 0.3 | 0.3 | 0.31 | 0.32 | 0.32 |
| 8 | mort | Residual | 0.46 | 0.49 | 0.52 | 0.57 | 0.6 |
| 9 | rec | quadrat | 0.18 | 0.15 | 0.15 | 0.11 | 0.1 |
| 10 | rec | quadrat:sp | 0.16 | 0.14 | 0.12 | 0.09 | 0.08 |
| 11 | rec | sp | 0.39 | 0.4 | 0.4 | 0.43 | 0.43 |
| 12 | rec | Residual | 0.27 | 0.32 | 0.34 | 0.37 | 0.39 |
Figure 1 complete
map
###### xy locations ######
forest <- read.table(here("data/plot_sites_information.csv"), sep=",",
header=T)
forest$long[forest$ID == "len"] <- forest$long[forest$ID == "len"] -2
fxy <- st_as_sf(forest[,-27], coords = c("long", "lat"), crs = 4326) #WSG84
fxy <- cbind(fxy, forest[,9:10])
fxy$lat = st_coordinates(fxy)[,2]
fxy$long = st_coordinates(fxy)[,1]
fxy$site[fxy$site == "Smithsonian Environmental Research Center"] <- "SERC"
fxy$site[fxy$site == "Smithsonian Conservation Biology Institute"] <- "SCBI"
world <- ne_countries(scale = "medium", returnclass = "sf") %>%
filter(geounit != "Antarctica") set.seed(5)
A <- ggplot(data = world) +
geom_sf(size=0.2) +
geom_sf(data=fxy, aes(fill=abs(lat), size=n_census_intervals),
shape=21, col="black") +
scale_size(name="# census \n intervals",range=c(2,10), breaks = c(1,2,3,5,7),
)+
scale_fill_gradientn(colours=pal) +
guides(fill = "none")+
geom_label_repel(data=fxy, aes(y=lat, x=long, label=site), size=5.5,
label.size = 0,box.padding = 0.7,
label.padding = 0.01) +
scale_x_continuous(expand = expansion(mult=c(-0.12, -0.15))) +
scale_y_continuous(expand = expansion(mult=-c(0.02, 0.15))) +
xlab("") + ylab("") +
theme_bw()+
theme(legend.position= c(0.05,0.18),
legend.text = element_text(size=12),
legend.title = element_text(size=15),
axis.text = element_blank(),
axis.ticks = element_blank(),
plot.margin = margin(-0.5, 0.1, -1, -0.3, "cm")) + labs(tag="(a)") +
geom_hline(yintercept=c(23.5,0,-23.5), linetype="dashed", col="gray") +#+
theme(plot.tag.position = c(0,1),
plot.tag = element_text(size=20, vjust = 1.2, hjust = 0, face="bold"))
B <- vpc.fig1 + ggtitle("") + labs(tag = "(b)") +
theme_cowplot()+
theme(plot.margin = margin(-0.5, 0.1, 0.2, 0, "cm"),
legend.position= "none",
#plot.tag.position = c(0,1),
# plot.tag = element_text(vjust = 3, hjust = -1.2),
text = element_text(size=20),
axis.text = element_text(size=16))
C <- tot.vpc.fig + labs(tag = "(c)") +
theme_cowplot()+
theme(plot.margin = margin(0.3, 0, 0.2, 0.3, "cm"),
legend.position= "none",
#plot.tag.position = c(0,1),
panel.spacing.x = unit(0.01, "cm"),
# plot.tag = element_text(vjust = 3, hjust = 0),
text = element_text(size=20),
axis.ticks = element_blank(),
axis.text.y = element_text(size=18),
axis.text.x = element_blank())
bot <- plot_spacer() + B + C + plot_layout(widths = c(0.00005,1,0.6))
#ggsave("figs/FIG_1BC.jpeg", height = 5.5, width=15)
A + wrap_plots(bot)+
plot_layout(ncol=1,widths = c(0.1,1), heights = c(1,0.6))ggsave(here("figures/FIG_1.jpeg"), height = 12, width=15)DIRICHLET regression
Including latitude for color
load(here("data", "plots_structure.Rdata"))
vital <- vital %>% left_join(plots.structure[,1:2], "fplot")Growth
predirig <- vital %>%
filter(vital == "grow", q_size == "5") %>%
select(fplot, time, term, VPC, richness.rarefaction, lat) %>%
ungroup() %>% select(-time) %>%
group_by(fplot, term) %>% summarise_all(mean)
diridatag <- predirig %>%
pivot_wider(names_from = term, values_from = VPC) %>% ungroup() %>%
mutate(log.rich = log(richness.rarefaction),
log.rich.o=log.rich) %>%
mutate_at(vars(log.rich), scale)
vpcg <- DR_data(diridatag[,c("quadrat", "sp", "quadrat:sp", "Residual")])
#plot(vpcg)mgrow <- DirichReg(vpcg~log.rich, data=diridatag, model="alternative", base=4)
summary(mgrow)## Call:
## DirichReg(formula = vpcg ~ log.rich, data = diridatag, model = "alternative",
## base = 4)
##
## Standardized Residuals:
## Min 1Q Median 3Q Max
## quadrat -0.8965 -0.5134 -0.2588 -0.0935 0.8402
## sp -1.7230 -1.0306 -0.0084 0.6703 3.3811
## quadrat:sp -1.5390 -0.6531 -0.2894 0.4352 3.2936
## Residual -2.1805 -0.8043 0.0819 0.7461 2.1985
##
## MEAN MODELS:
## ------------------------------------------------------------------
## Coefficients for variable no. 1: quadrat
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.5532 0.1760 -14.50 <2e-16 ***
## log.rich 0.4056 0.1711 2.37 0.0178 *
## ------------------------------------------------------------------
## Coefficients for variable no. 2: sp
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.63578 0.08934 -7.116 1.11e-12 ***
## log.rich -0.28094 0.09278 -3.028 0.00246 **
## ------------------------------------------------------------------
## Coefficients for variable no. 3: quadrat:sp
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.37320 0.11387 -12.059 <2e-16 ***
## log.rich -0.04767 0.11216 -0.425 0.671
## ------------------------------------------------------------------
## Coefficients for variable no. 4: Residual
## - variable omitted (reference category) -
## ------------------------------------------------------------------
##
## PRECISION MODEL:
## ------------------------------------------------------------------
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.4508 0.1791 19.27 <2e-16 ***
## ------------------------------------------------------------------
## Significance codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-likelihood: 102.9 on 7 df (45 BFGS + 1 NR Iterations)
## AIC: -191.7, BIC: -184.4
## Number of Observations: 21
## Links: Logit (Means) and Log (Precision)
## Parametrization: alternative
Predictions
newdata <- expand.grid(log.rich = seq(min(diridatag$log.rich),
max(diridatag$log.rich), length.out=10))
pred <- predict(mgrow, newdata = newdata, se=T)
colnames(pred) <- colnames(vpcg)
newdata$log.rich.o <- newdata$log.rich*sd(diridatag$log.rich.o) +
mean(diridatag$log.rich.o)
newdata$rich.o <- exp(newdata$log.rich.o)
newpredg <- cbind(newdata,pred) %>% pivot_longer(`quadrat`:`Residual`, names_to="term", values_to="pred")
newpredg %>% mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp",
"Residual")) %>%
ggplot(aes(x=rich.o, y=pred, col=term)) +
geom_line() + facet_grid(~term) +
geom_point(data=predirig, aes(x=richness.rarefaction, y=VPC, col=term)) +
scale_x_log10() +
ggtitle("Mod log.rich")
Residuals
resid <- residuals(mgrow, type = "standardized")
resid <- data.frame(resid[,1:4])
resid$log.rich <- diridatag$log.rich
resid <- resid %>% pivot_longer(1:4, names_to="term", values_to="resid")
ggplot(resid, aes(x=log.rich, y=resid)) + geom_point() +
facet_grid(~term)Mortality
predirim <- vital %>%
filter(vital == "mort", q_size == "5") %>%
select(fplot, time, term, VPC,richness.rarefaction, lat) %>%
ungroup() %>% select(-time) %>%
group_by(fplot, term) %>% summarise_all(mean)
diridatam <- predirim %>%
pivot_wider(names_from = term, values_from = VPC) %>% ungroup() %>%
mutate(log.rich = log(richness.rarefaction),
log.rich.o=log.rich) %>%
mutate_at(vars(log.rich), scale)
vpcm <- DR_data(diridatam[,c("quadrat", "sp", "quadrat:sp", "Residual")])mmort <- DirichReg(vpcm~log.rich,data=diridatam, model="alternative", base=4)
summary(mmort)## Call:
## DirichReg(formula = vpcm ~ log.rich, data = diridatam, model = "alternative",
## base = 4)
##
## Standardized Residuals:
## Min 1Q Median 3Q Max
## quadrat -1.2381 -0.7293 -0.0222 0.2704 3.7743
## sp -2.3194 -0.4716 0.3618 0.7887 3.1973
## quadrat:sp -1.4603 -0.5624 -0.0776 0.1208 2.1838
## Residual -2.0367 -0.5814 0.0805 0.5724 1.9470
##
## MEAN MODELS:
## ------------------------------------------------------------------
## Coefficients for variable no. 1: quadrat
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.6337 0.1444 -11.311 <2e-16 ***
## log.rich 0.3434 0.1564 2.196 0.0281 *
## ------------------------------------------------------------------
## Coefficients for variable no. 2: sp
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.56723 0.09903 -5.728 1.02e-08 ***
## log.rich -0.11602 0.09819 -1.182 0.237
## ------------------------------------------------------------------
## Coefficients for variable no. 3: quadrat:sp
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.2057 0.1224 -9.848 <2e-16 ***
## log.rich -0.2129 0.1315 -1.618 0.106
## ------------------------------------------------------------------
## Coefficients for variable no. 4: Residual
## - variable omitted (reference category) -
## ------------------------------------------------------------------
##
## PRECISION MODEL:
## ------------------------------------------------------------------
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.2755 0.1752 18.69 <2e-16 ***
## ------------------------------------------------------------------
## Significance codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-likelihood: 85.36 on 7 df (44 BFGS + 1 NR Iterations)
## AIC: -156.7, BIC: -149.4
## Number of Observations: 21
## Links: Logit (Means) and Log (Precision)
## Parametrization: alternative
newdata <- expand.grid(log.rich = seq(min(diridatam$log.rich),
max(diridatam$log.rich), length.out=10))
pred <- predict(mmort, newdata = newdata, se=T)
confint(mmort)##
## 95% Confidence Intervals (original form)
##
## - Beta-Parameters:
## Variable: quadrat
## 2.5% Est. 97.5%
## (Intercept) -1.917 -1.634 -1.35
## log.rich 0.037 0.343 0.65
##
## Variable: sp
## 2.5% Est. 97.5%
## (Intercept) -0.761 -0.567 -0.373
## log.rich -0.308 -0.116 0.076
##
## Variable: quadrat:sp
## 2.5% Est. 97.5%
## (Intercept) -1.446 -1.206 -0.966
## log.rich -0.471 -0.213 0.045
##
## Variable: Residual
## variable omitted
##
## - Gamma-Parameters
## 2.5% Est. 97.5%
## (Intercept) 2.93 3.27 3.62
colnames(pred) <- colnames(vpcm)
newdata$log.rich.o <- newdata$log.rich*sd(diridatam$log.rich.o) +
mean(diridatam$log.rich.o)
newdata$rich.o <- exp(newdata$log.rich.o)
newpredm <- cbind(newdata,pred) %>% pivot_longer(`quadrat`:`Residual`, names_to="term", values_to="pred")
newpredm %>% mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp",
"Residual")) %>%
ggplot(aes(x=rich.o, y=pred, col=term)) +
geom_line() + facet_grid(~term) +
geom_point(data=predirim, aes(x=richness.rarefaction, y=VPC, col=term)) +
scale_x_log10() +
ggtitle("Mod log.rich")Residuals
resid <- residuals(mmort, type = "standardized")
resid <- data.frame(resid[,1:4])
resid$log.rich <- diridatam$log.rich
resid <- resid %>% pivot_longer(1:4, names_to="term", values_to="resid")
ggplot(resid, aes(x=log.rich, y=resid)) + geom_point() +
facet_grid(~term)Recruitment
predirir <- vital %>%
filter(vital == "rec", q_size == "5") %>%
select(fplot, time, term, VPC,richness.rarefaction, lat) %>%
ungroup() %>% select(-time) %>%
group_by(fplot, term) %>% summarise_all(mean)
diridatar<- predirir %>%
pivot_wider(names_from = term, values_from = VPC) %>% ungroup() %>%
mutate(log.rich = log(richness.rarefaction),
log.rich.o=log.rich) %>%
mutate_at(vars(log.rich), scale)
vpcr <- DR_data(diridatar[,c("quadrat", "sp", "quadrat:sp", "Residual")])mrec <- DirichReg(vpcr~log.rich,data=diridatar, model="alternative", base=4)
summary(mrec)## Call:
## DirichReg(formula = vpcr ~ log.rich, data = diridatar, model = "alternative",
## base = 4)
##
## Standardized Residuals:
## Min 1Q Median 3Q Max
## quadrat -1.5135 -0.6182 -0.1633 0.6980 2.0837
## sp -2.8000 -0.5444 0.0109 0.5788 1.9356
## quadrat:sp -1.5556 -0.6095 0.0688 0.4667 2.5132
## Residual -1.6572 -0.6812 0.0204 0.3700 2.4697
##
## MEAN MODELS:
## ------------------------------------------------------------------
## Coefficients for variable no. 1: quadrat
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.58720 0.12451 -4.716 2.4e-06 ***
## log.rich 0.03787 0.13509 0.280 0.779
## ------------------------------------------------------------------
## Coefficients for variable no. 2: sp
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.08644 0.10441 0.828 0.408
## log.rich -0.50953 0.10712 -4.757 1.97e-06 ***
## ------------------------------------------------------------------
## Coefficients for variable no. 3: quadrat:sp
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.6866 0.1273 -5.396 6.83e-08 ***
## log.rich -0.1179 0.1362 -0.865 0.387
## ------------------------------------------------------------------
## Coefficients for variable no. 4: Residual
## - variable omitted (reference category) -
## ------------------------------------------------------------------
##
## PRECISION MODEL:
## ------------------------------------------------------------------
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.3645 0.1829 18.4 <2e-16 ***
## ------------------------------------------------------------------
## Significance codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-likelihood: 73.9 on 7 df (39 BFGS + 1 NR Iterations)
## AIC: -133.8, BIC: -127.2
## Number of Observations: 19
## Links: Logit (Means) and Log (Precision)
## Parametrization: alternative
Prediction
newdata <- expand.grid(log.rich = seq(min(diridatar$log.rich),
max(diridatar$log.rich), length.out=10))
pred <- predict(mrec, newdata = newdata, se=T)
colnames(pred) <- colnames(vpcr)
newdata$log.rich.o <- newdata$log.rich*sd(diridatar$log.rich.o) +
mean(diridatar$log.rich.o)
newdata$rich.o <- exp(newdata$log.rich.o)
newpredr <- cbind(newdata,pred) %>% pivot_longer(`quadrat`:`Residual`, names_to="term", values_to="pred")
newpredr %>% mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp",
"Residual")) %>%
ggplot(aes(x=rich.o, y=pred, col=term)) +
geom_line() + facet_grid(~term) +
geom_point(data=predirir, aes(x=richness.rarefaction, y=VPC, col=term)) +
scale_x_log10() +
ggtitle("Mod log.rich")Residuals
resid <- residuals(mrec, type = "standardized")
resid <- data.frame(resid[,1:4])
resid$log.rich <- diridatar$log.rich
resid <- resid %>% pivot_longer(1:4, names_to="term", values_to="resid")
ggplot(resid, aes(x=log.rich, y=resid)) + geom_point() +
facet_grid(~term)Figure dirichlet
Prediction intervals calculated in another code following the paper:
Douma, J.C. & Weedon, J.T. (2019) Analysing continuous proportions in ecology and evolution: A practical introduction to beta and Dirichlet regression. Methods in Ecology and Evolution, 10, 1412–1430. https://doi.org/10.1111/2041-210X.13234.
source(here("scripts", "prediction_intervals_dirichlet.R"), local=T)predis <- bind_rows(list(grow=newpredg, mort=newpredm, rec=newpredr),
.id="vital") %>%
mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp","Residual"))
pontos <- bind_rows(list(grow=predirig, mort= predirim, rec=predirir), .id="vital")%>%
mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp", "Residual"))
quants <- bind_rows(list(grow=quantg,mort=quantm, rec=quantr), .id="vital")%>%
mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp", "Residual"))# pvalues
res <- data.frame(vital = c("grow", "mort", "rec"),
P = NA,
term = "Residual")
test <- summary(mgrow)
vals <- as.data.frame(test$coef.mat)[1:6,]
vals$coef <- names(test$coefficients)[1:6]
valsg <- vals[ grep("log.rich", vals$coef),] %>%
mutate(term=c("quadrat", "sp", "quadrat:sp"))
test <- summary(mmort)
vals <- as.data.frame(test$coef.mat)[1:6,]
vals$coef <- names(test$coefficients)[1:6]
valsm <- vals[ grep("log.rich", vals$coef),] %>%
mutate(term=c("quadrat", "sp", "quadrat:sp"))
test <- summary(mrec)
vals <- as.data.frame(test$coef.mat)[1:6,]
vals$coef <- names(test$coefficients)[1:6]
valsr <- vals[ grep("log.rich", vals$coef),] %>%
mutate(term=c("quadrat", "sp", "quadrat:sp"))
pvals <- bind_rows(list(grow=valsg,mort=valsm,rec=valsr), .id="vital") %>%
dplyr::select(vital,`Pr(>|z|)`, term) %>% rename(P = `Pr(>|z|)`)
pvals <- rbind(pvals,res) %>%
mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp", "Residual"))
pvals$x <- 300
pvals$y <- 0.70
pvals$sig <- "ns"
pvals$sig[which(pvals$P < 0.016667)] <- paste0("p = ", round(pvals$P[which(pvals$P < 0.016667)],3))
pvals$sig[pvals$term=="Residual"] <- ""
pvals$sig[pvals$sig=="p = 0"] <- "p < 0.001"cfs <- bind_rows(list(grow=valsg,mort=valsm,rec=valsr), .id="vital") %>%
select(vital, coef,term, Estimate, `Std. Error`,`Pr(>|z|)`)
ggplot(cfs, aes(x=Estimate, y=term, col=term)) +
facet_grid(~vital) +
geom_point() +
geom_linerange(aes(xmin=Estimate-`Std. Error`, xmax=Estimate+`Std. Error`))+
geom_vline(xintercept = 0, linetype="dashed") +
ggtitle("Log.Rich Estimates + Std.Error")library(wesanderson) #better colour palette
# Gradient color para latitutde
pal <- wes_palette("Zissou1",21, type = "continuous")[21:1] # azul frio-temperadofdiri_lat <-ggplot(predis, aes(x=rich.o, y=pred)) +
geom_line(size=1) +
facet_grid(vital~term, labeller=grpvit) +
geom_smooth(data=quants,aes(x=rich.o, y=lower), se=F, size=0.1)+
geom_smooth(data=quants,aes(x=rich.o, y=upper), se=F, size=0.1) +
geom_ribbon(data=quants,aes(x=rich.o, ymin=lower, ymax=upper,
y=mean), alpha=0.05, fill="blue",
size=0)+
geom_point(data=pontos, aes(x=richness.rarefaction, y=VPC, fill=abs(lat)),
col="black", size=2.5, pch=21)+
scale_fill_gradientn(colors=pal) +
scale_y_continuous(minor_breaks = NULL) +
scale_x_log10() +
geom_text(data=pvals, aes(x=x,y=y, label=sig), size=3, hjust=1,vjust=0,
col="black")+
xlab("Species richness")+
ylab("VPC") +
theme_cowplot() +
theme(panel.background = element_rect(colour="lightgray"),
legend.position = "none")
fdiri_latpng(here("figures", "FIG_4_dirichlet_regressions.png"), height = 600, width=800, res=100)
fdiri_lat
dev.off()## quartz_off_screen
## 2